home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 2.iso / BARNET / COMPILER / SATHER / !Sather / Library / Math / sa / svd < prev    next >
Text File  |  1996-07-25  |  12KB  |  341 lines
  1. ---------------------------> Sather 1.1 source file <--------------------------
  2. -- COPYRIGHT NOTICE: This code is provided WITHOUT ANY WARRANTY
  3. -- and is subject to the terms of the SATHER LIBRARY GENERAL PUBLIC
  4. -- LICENSE contained in the file: Sather/Doc/License of the
  5. -- Sather distribution. The license is also available from ICSI,
  6. -- 1947 Center St., Suite 600, Berkeley CA 94704, USA.
  7. -------------------------------------------------------------------
  8. -- This class had calls to a c version that seems to no longer exist!
  9. -- I commented that stuff out and added the test to test-oblig.sa
  10. -- (gomes Nov21/95)
  11. -- *************************************************************
  12. -- Converted to use the new matrix and vector classes July 96
  13. -- The test class for SVD seems to generate results that are slightly
  14. -- off from the expected results in 0.2. 
  15. -- Bo Hagstedt<bo.hagstedt@mailbox.swipnet.se> looked at the code
  16. -- and found that it was good (and fixed several errors).
  17. -- The IF statments protect accesses outside the matrices when
  18. -- the loop index i points to the last row column as the sub loop will never 
  19. -- start when the start index is greater han the stop indexand a sub loop uses
  20. -- i+1 as its index.  
  21. -- Thanks also to Alex Cozzi who looked at this code earlier.
  22. -- gomes@icsi.berkeley.edu
  23. -------------------------------------------------------------------
  24. class NR_SVD is
  25.    -- Singular value decomposition derived from Numerical Recipes.
  26.    -- This algorithm is known to have a problem with an obscure case (but
  27.    -- works almost all of the time). Better algorithms are being worked on.
  28.    
  29.    pythag(a,b:FLT):FLT is
  30.       -- Square root of a^2+b^2 without destructive overflow or underflow.
  31.       res:FLT;
  32.       at:FLT:=a.abs;  bt:FLT:=b.abs; ct:FLT;
  33.       if at>bt then ct:=bt/at; res:=at*(1.0+ct*ct).sqrt;
  34.       elsif bt/=0.0 then ct:=at/bt; res:=bt*(1.0+ct*ct).sqrt;
  35.       end; -- if
  36.       return res;
  37.    end; -- pythag
  38.    
  39.    svd(a:MAT, w:VEC, v:MAT)
  40.       pre a.nr>=a.nc and w.dim=a.nc and v.nr=a.nc and v.nc=a.nc    is
  41.      -- Computes the singular value decomposition of a. When finished,
  42.      -- a w v^T equals the old a. Must have a.nr>=a.nc.
  43.      -- Based on Numerical Recipes in C, p. 68.
  44.       m:INT:=a.nr; n:INT:=a.nc;
  45.       -- Householder reduction to bidiagonal form (should be separate).
  46.       rv1:VEC:=rv1.create(n); 
  47.       g,scale,anorm,c,f,h,s,x,y,z:FLT; 
  48.       i,its,j,jj,k,l,nm:INT;  flag:BOOL;
  49.       -- Householder reduction to bidiagonal form
  50.       i:=0; loop until!(i>=n);
  51.      l:=i+1;  rv1[i]:=scale*g; g:=0.0; s:=0.0; scale:=0.0;
  52.      if i<m then 
  53.         k:=i; loop until!(k>=m); scale:=scale+a[k,i].abs; k:=k+1 end; 
  54.         if scale/=0.0 then
  55.            k:=i; loop until!(k>=m); 
  56.           a[k,i]:=a[k,i]/scale; s:=s+a[k,i]*a[k,i]; k:=k+1 
  57.            end; -- loop
  58.            f:=a[i,i]; 
  59.            if f<0.0 then g:=s.sqrt else g:=-s.sqrt end;
  60.            h:=f*g-s; a[i,i]:=f-g;
  61.            j:=l; loop until!(j>=n);
  62.           s:=0.0; 
  63.           k:=i; loop until!(k>=m); s:=s+a[k,i]*a[k,j]; k:=k+1 end;
  64.           f:=s/h;
  65.           k:=i; loop until!(k>=m); a[k,j]:=a[k,j]+f*a[k,i];k:=k+1 end;
  66.           j:=j+1
  67.            end; -- loop
  68.            k:=i; loop until!(k>=m); a[k,i]:=a[k,i]*scale; k:=k+1 end;
  69.         end; -- if
  70.      end; -- if
  71.      w[i]:=scale*g; g:=0.0; s:=0.0; scale:=0.0;
  72.      if i<m and i/=n-1 then
  73.         k:=l; loop until!(k>=n); scale:=scale+a[i,k].abs; k:=k+1 end;
  74.         if scale/=0.0 then
  75.            k:=l; loop until!(k>=n);
  76.           a[i,k]:=a[i,k]/scale; s:=s+a[i,k]*a[i,k]; k:=k+1
  77.            end; -- loop
  78.            f:=a[i,l];
  79.            if f<0.0 then g:=s.sqrt else g:=-s.sqrt end;
  80.            h:=f*g-s; a[i,l]:=f-g;
  81.            k:=l; loop until!(k>=n); 
  82.           rv1[k]:=a[i,k]/h;
  83.           k:=k+1 
  84.            end;
  85.            j:=l; loop until!(j>=m); 
  86.           s:=0.0;
  87.           k:=l; loop until!(k>=n); s:=s+a[j,k]*a[i,k]; k:=k+1 end;
  88.           k:=l; loop until!(k>=n); a[j,k]:=a[j,k]+s*rv1[k];k:=k+1 end;
  89.           j:=j+1
  90.            end; -- loop
  91.            k:=l; loop until!(k>=n); a[i,k]:=a[i,k]*scale; k:=k+1 end;
  92.         end; -- if
  93.      end; -- if
  94.      anorm:=anorm.max(w[i].abs+rv1[i].abs);
  95.      i:=i+1;
  96.       end; -- loop
  97.       -- Accumulation of right-hand transformations.
  98.       i:=n-1; loop until!(i<0); 
  99.      if i<n-1 then
  100.         if g/=0.0 then
  101.            j:=l; loop until!(j>=n); 
  102.           -- double division to avoid possible underflow
  103.           v[j,i]:=(a[i,j]/a[i,l])/g; 
  104.           j:=j+1 
  105.            end;
  106.            j:=l; loop until!(j>=n);
  107.           s:=0.0;
  108.           k:=l; loop until!(k>=n); s:=s+a[i,k]*v[k,j]; k:=k+1 end;
  109.           k:=l; loop until!(k>=n); v[k,j]:=v[k,j]+s*v[k,i]; k:=k+1 end;
  110.           j:=j+1
  111.            end; -- loop
  112.         end; -- if
  113.         j:=l; loop until!(j>=n); v[i,j]:=0.0; v[j,i]:=0.0; j:=j+1 end;
  114.      end; -- if
  115.      v[i,i]:=1.0; g:=rv1[i]; l:=i;
  116.      i:=i-1;
  117.       end; -- loop
  118.       -- Accumulation of left-hand transformations.
  119.       -- the new code is different : i=IMIN(m,n)
  120.       --      i:=n-1; loop until!(i<0); 
  121.       i:=n.min(m)-1; loop until!(i<0); 
  122.      l:=i+1; g:=w[i];
  123.      -- Redundant if removed by BH
  124.      j:=l; loop until!(j>=n); a[i,j]:=0.0; j:=j+1 end;
  125.      if g/=0.0 then
  126.         g:=1.0/g;
  127.         
  128.         -- if i/=n-1 then BH
  129.         j:=l; loop until!(j>=n);
  130.            s:=0.0;
  131.            k:=l; loop until!(k>=m); s:=s+a[k,i]*a[k,j]; k:=k+1 end;
  132.            f:=(s/a[i,i])*g;
  133.            k:=i; loop until!(k>=m); a[k,j]:=a[k,j]+f*a[k,i]; k:=k+1 end;
  134.            j:=j+1
  135.         end; -- loop
  136.         -- end; -- if
  137.         j:=i; loop until!(j>=m); a[j,i]:=a[j,i]*g; j:=j+1 end;
  138.      else
  139.         j:=i; loop until!(j>=m); a[j,i]:=0.0; j:=j+1 end;
  140.      end; -- if
  141.      a[i,i]:=a[i,i]+1.0;
  142.      i:=i-1;
  143.       end; -- loop
  144.       -- Diagonalization of the bidiagonal form.
  145.       k:=n-1; loop until!(k<0);    -- loop over singular values
  146.      --its:=0; loop until!(its>=30); -- loop over allowed iterations
  147.      -- Should read as follows in order to get the "its" error message. /BH
  148.      its:=1; loop until!(its>30); -- loop over allowed iterations
  149.         flag:=true; 
  150.         l:=k; loop until!(l<0); -- test for splitting
  151.            nm:=l-1;        -- rv[1] is always zero
  152.            if rv1[l].abs.flt+anorm=anorm then flag:=false; break! end;
  153.            if w[nm].abs.flt+anorm=anorm then break! end;
  154.            l:=l-1
  155.         end; -- loop
  156.         if flag then
  157.            c:=0.0; s:=1.0; 
  158.            i:=l; loop until!(i>k);
  159.           f:=s*rv1[i];
  160.           -- Corrected the typo, rv1[1] => rv1[i]. /BH
  161.           rv1[i]:=c*rv1[i];
  162.           if f.abs.flt+anorm/=anorm then
  163.              g:=w[i]; h:=pythag(f,g); w[i]:=h; h:=1.0/h;
  164.              c:=g*h; s:=-f*h;
  165.              j:=0; loop until!(j>=m); 
  166.             y:=a[j,nm]; z:=a[j,i]; 
  167.             a[j,nm]:=y*c+z*s; a[j,i]:=z*c-y*s;
  168.             j:=j+1
  169.              end; -- loop
  170.           end; -- if
  171.           i:=i+1
  172.            end; -- loop
  173.         end; -- if
  174.         z:=w[k];
  175.         if l=k then        -- convergence
  176.            if z<0.0 then    -- singular value is made non-negative
  177.           w[k]:=-z;
  178.           j:=0; loop until!(j>=n); v[j,k]:=-v[j,k]; j:=j+1 end;
  179.            end; -- if
  180.            break!;
  181.         end; -- if
  182.         if its=30 then #ERR+"SVD: 30 iterations, no convergence\n"; end;
  183.         -- shift from bottom 2 by 2 minor
  184.         x:=w[l]; nm:=k-1; y:=w[nm]; g:=rv1[nm]; h:=rv1[k];
  185.         f:=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); g:=pythag(f,1.0);
  186.         if f>=0.0 then f:=((x-z)*(x+z)+h*((y/(f+g.abs))-h))/x
  187.         else f:=((x-z)*(x+z)+h*((y/(f-g.abs))-h))/x end;
  188.         -- Next QR transformation
  189.         c:=1.0; s:=1.0; 
  190.         j:=l; loop until!(j>=nm+1);
  191.            i:=j+1; g:=rv1[i]; y:=w[i]; h:=s*g; g:=c*g; z:=pythag(f,h);
  192.            rv1[j]:=z; c:=f/z; s:=h/z; f:=x*c+g*s; g:=g*c-x*s;
  193.            h:=y*s; y:=y*c; 
  194.            jj:=0; loop until!(jj>=n);
  195.           x:=v[jj,j]; z:=v[jj,i]; v[jj,j]:=x*c+z*s; v[jj,i]:=z*c-x*s;
  196.           jj:=jj+1
  197.            end; -- loop
  198.            z:=pythag(f,h); w[j]:=z;
  199.            if z/=0.0 then z:=1.0/z; c:=f*z; s:=h*z end;
  200.            f:=(c*g)+(s*y); x:=(c*y)-(s*g);
  201.            jj:=0; loop until!(jj>=m); 
  202.           y:=a[jj,j]; z:=a[jj,i]; a[jj,j]:=y*c+z*s; a[jj,i]:=z*c-y*s;
  203.           jj:=jj+1;
  204.            end; -- loop
  205.            j:=j+1
  206.         end; -- loop
  207.         rv1[l]:=0.0; rv1[k]:=f; w[k]:=x;
  208.         its:=its+1;
  209.      end; -- loop
  210.      k:=k-1
  211.       end; -- loop
  212.    end; -- svd
  213.  
  214.    
  215. end; -- class NR_SVD
  216.  
  217. -------------------------------------------------------------------
  218.  
  219. class TEST_SVD is
  220.    include TEST;
  221.  
  222.    main is
  223.       class_name("SVD");
  224.       a:MAT:=MAT::create(3,4);
  225.       a.inplace_uniform_random;
  226.       svd_matrix_test(1, a);
  227.       a:=MAT::create(3,4);
  228.       a.inplace_uniform_random;
  229.       svd_matrix_test(2, a);
  230.       a:=MAT::create(4,3); a.inplace_uniform_random;
  231.       svd_matrix_test(3, a);
  232.       a:=MAT::create(5, 5); a.inplace_uniform_random;
  233.       svd_matrix_test(4, a);
  234.       a:=MAT::create(||0.0,1.0,0.0|,|0.0,1.0,1.0|,|0.0,0.0,0.0||);
  235.       svd_matrix_test(5, a);    -- the bad example for Numerical Recipes
  236.  
  237.       svd_back_sub_test(1,
  238.         MAT::create(||1.0,0.0,0.0|,|0.0,2.0,0.0|,|0.0,0.0,3.0||),
  239.         VEC::create(|1.0,0.0,0.0|));
  240.       svd_back_sub_test(2,
  241.         MAT::create(||1.0,0.0,0.0|,|0.0,2.0,0.0|,|0.0,0.0,3.0||),
  242.         VEC::create(|1.0,1.0,1.0|));
  243.       svd_back_sub_test(3,
  244.         MAT::create(||1.0,2.0,3.0|,|2.0,2.0,3.0|,|3.0,2.0,3.0||),
  245.         VEC::create(|1.0,1.0,1.0|));
  246.       svd_back_sub_test(4,
  247.         MAT::create(||1.0,2.0,3.0|,|2.0,2.0,3.0|,|3.0,2.0,3.0||),
  248.         VEC::create(|1.0,2.0,3.0|));
  249.       svd_back_sub_test(5,
  250.         MAT::create(||1.0,2.0,3.0,4.0,5.0|,|2.0,3.0,4.0,5.0,1.0|
  251.         ,|3.0,4.0,5.0,1.0,2.0|,|4.0,5.0,1.0,2.0,3.0|,|5.0,1.0,2.0,3.0,4.0||),
  252.         VEC::create(|1.0,1.0,1.0,1.0,1.0|));
  253.       svd_back_sub_test(6,
  254.         MAT::create(||1.0,2.0,3.0,4.0,5.0|,|2.0,3.0,4.0,5.0,1.0|
  255.         ,|3.0,4.0,5.0,1.0,2.0|,|4.0,5.0,1.0,2.0,3.0|,|5.0,1.0,2.0,3.0,4.0||),
  256.         VEC::create(|1.0,2.0,3.0,4.0,5.0|));
  257.       svd_back_sub_test(7,
  258.         MAT::create(||1.4, 2.1, 2.1, 7.4, 9.6|
  259.         , |1.6, 1.5, 1.1, 0.7, 5.0|, |3.8, 8.0, 9.6, 5.4, 8.8|
  260.         , |4.6, 8.2, 8.4, 0.4, 8.0|, |2.6, 2.9, 0.1, 9.6, 7.7||),
  261.         VEC::create(|1.1, 1.6, 4.7, 9.1, 0.1|));
  262.  
  263.       lfm:MAT:=MAT::create(||1.0,2.0|,|3.0,4.0|,|5.0,6.0||);
  264.       fl_lfi:FLIST{VEC}:=FLIST{VEC}::create;
  265.       fl_lfi:=fl_lfi.push(VEC::create(|1.0,0.0|));
  266.       fl_lfi:=fl_lfi.push(VEC::create(|1.0,1.0|));
  267.       fl_lfi:=fl_lfi.push(VEC::create(|0.0,1.0|));
  268.       fl_lfi:=fl_lfi.push(VEC::create(|2.0,1.0|));
  269.       lfi:ARRAY{VEC} := fl_lfi.array;
  270.       
  271.       fl_lfo:FLIST{VEC}:=FLIST{VEC}::create; 
  272.       loop
  273.      fl_lfo:=fl_lfo.push(lfm.times_vec(lfi.elt!));
  274.       end;
  275.       lfo: ARRAY{VEC} := fl_lfo.array;
  276.       loop
  277.      #OUT +"lfo: "+lfo.elt!.str+"\n";
  278.       end;
  279.       
  280. --d      test("to_linear_fit_of",
  281. --d        lfm.copy.inplace_zero.to_linear_fit_of(lfi,lfo).str, lfm.str);
  282.       fl_lfoa:FLIST{VEC}:=FLIST{VEC}::create;
  283.       loop
  284.      fl_lfoa:=fl_lfoa.push(lfo.elt!+(VEC::create(|1.0,2.0,3.0|)));
  285.       end;
  286.       lfoa:ARRAY{VEC} := fl_lfoa.array;
  287.       
  288.       lfma:MAT:=MAT::create
  289.         (||1.0,2.0,1.0|,|3.0,4.0,2.0|,|5.0,6.0,3.0||);
  290. --d      test("to_affine_fit_of",
  291. --d        lfma.copy.inplace_zero.inplace_affine_fit_of(lfi,lfoa).str, lfma.str);
  292. --d      wt:FLIST{FLT}:=FLIST{FLT}::create;
  293. --d      wt:=wt.push(1.0); wt:=wt.push(1.0); wt:=wt.push(1.0); wt:=wt.push(1.0);
  294. --d      lfmcopy ::= lfm.copy;
  295. --d      lfmcopy.inplace_zero;
  296. --d      lfmcopy.inplace_weighted_linear_fit_of(lfi,lfo,wt);
  297. --d      test("to_weighted_linear_fit_of",lfmcopy.str,lfm.str);
  298. --d      test("to_weighted_affine_fit_of",
  299. --d        lfma.copy.inplace_zero.inplace_weighted_affine_fit_of(lfi,lfoa,wt).str
  300. --d        , lfma.str);
  301.       finish;
  302.    end;
  303.  
  304.    svd_matrix_test(n:INT, a:MAT) is
  305.       -- Test the singular value decomposition of `a' labelling the test
  306.       -- by `n'.
  307.       u:MAT:=u.create(a.nr.max(a.nc),a.nc);
  308.       v:MAT:=v.create(a.nc, a.nc); w:VEC:=w.create(a.nc);
  309.       a.svd_in(u,w,v);
  310.       u.inplace_times_diagonal(w);
  311.       v.inplace_trans;
  312.       tmp:MAT:=u.copy;
  313.       tmp.inplace_arg_times_arg(u,v);
  314.       acpy::= a.copy;
  315.       acpy.inplace_zero;
  316.       acpy.inplace_portion_of_arg(tmp);
  317.       cmp:MAT:=acpy;
  318.       unchecked_test("svd_in "+n, cmp.str, a.str);
  319.    end; -- svd_test_matrix
  320.  
  321.    svd_back_sub_test(n:INT, a:MAT, b:VEC) is
  322.       -- Test `svd_back_sub' on the problem `a.x=b', label with `n'.
  323.       size:INT:=a.nr;
  324.       u:MAT:=u.create(size,size); v:MAT:=v.create(size,size);
  325.       w:VEC:=w.create(size); a.svd_in(u,w,v);
  326.       wmax:FLT:=w.max_value; wmin:FLT:=wmax*(0.000001);
  327.       k:INT:=0; loop until!(k=size);
  328.      if w[k]<=wmin then w[k]:=0.0 end; k:=k+1 
  329.       end; -- loop
  330.       x:VEC:=VEC::create(size);
  331.       a.svd_back_sub(u,w,v,b,x);
  332.       c:VEC:=a.times_vec(x);
  333.       unchecked_test("svd_back_sub "+n, c.str, b.str);
  334.    end; -- svd_back_sub_test
  335.  
  336.  
  337.  
  338. end; -- class TEST_SVD
  339.  
  340. -------------------------------------------------------------------
  341.